!> two-dimensional matrix
!
!  This module contains the data structure and procedures necessary
!  for creating and handling the coefficient matrix K in Kx=b.  The
!  matrix is a five band matrix and is generated first into a
!  diagonal structure and is then stored in the compressed row
!  format.
!
!  J. Roberts, 6/6/10
module class_2d_matx

    use matx_const
    use class_input, only: indat

    implicit none

    contains

!************************************************************
!> buildme takes the results of the coefficient generator and
!  builds the necessary linear system components
!************************************************************
    subroutine buildk(K,inp,rhs,dat)

        implicit none
#include "finclude/petsc.h90"
        !> array of matrices for multigroup diffusion operator
        Mat, dimension(:), pointer   :: K
        type(indat), intent(inout)   :: inp
        type(rhsdat), intent(inout)  :: rhs
        type(coefdat), intent(inout) :: dat
        integer                      :: g, n
        double precision             :: diag(inp%nxfm*inp%nyfm,5,inp%numg)
        PetscInt                     :: ierr
        ! matrix dimension
        n = inp%nxfm*inp%nyfm
        ! allocate K
        allocate( K(inp%numg) )
        ! now initialize the individual componenents
        do g = 1, inp%numg
            call MatCreate(PETSC_COMM_WORLD,K(g),ierr)
            call MatSetType(K(g),MATAIJ,ierr)
            call MatSetSizes(K(g),         & ! Mat
                             PETSC_DECIDE, & ! number of local rows
                             PETSC_DECIDE, & ! number of local col
                             n,            & ! number of global rows
                             n,            & ! number of global col
                             ierr          ) ! error
            call MatSetFromOptions(K(g),ierr)
        end do

        ! build the matrix
        call setk(K,inp,rhs,dat)
        write(0,*) " all done setting K "

    end subroutine buildk

    subroutine DestroyK(inp,K)
        implicit none
#include "finclude/petsc.h90"
        type(indat), intent(in)    :: inp
        Mat, dimension(:), pointer :: K
        PetscInt                   :: ierr, g
        do g = 1, inp%numg
            call MatDestroy(K(g),ierr)
        end do
        deallocate(K)
    end subroutine DestroyK

end module class_2d_matx
